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1.0  ACTIVITY  SUMMARY 


This  is  our  Phase  I  final  report  covering  the  period  from  04/08/05  to  01/07/06.  During  this  Phase 
I  period,  the  following  research  and  development  (R  &  D)  activities  have  taken  place: 

1.1  Revision  of  Phase  I  Statement  of  Works 

Based  on  the  suggestions  given  by  our  technical  monitors  during  the  kick-off  meeting,  we  have 
revised  our  Phase  I  statement  of  works.  The  major  tasks  of  our  Phase  I  contract  will  be  revised  to 
the  following  R&D: 

•  Improve  the  capability  of  the  current  Navier-Stokes  solvers  for  prediction  of  low-Reynolds- 
number  (low-Re)  airfoil  aerodynamics  by  adding  the  transition  prediction  capability  to  the 
Navier-Stokes  solvers. 

•  Validate  the  improved  Navier-Stokes  solver  with  Air  Force  Research  Laboratory’s  (AFRL) 
low-Re  SD7003  airfoil  case. 

1.2  Investigation  of  Low-Re  Version  of  Wilcox’s  k-<v  Model 

Based  on  our  previous  experience,  we  first  implement  and  test  low-Re  version  of  Wilcox’s  k-co 
model  for  prediction  of  laminar  separation  bubble.  Some  preliminary  results  are  presented  in 
Section  3.0. 

1.3  Reynolds  Averaged  Navier  Stokes  (RANS)  Simulation  of  Low-Re  Airfoil 
Aerodynamics 

After  the  above  warm-up  exercise,  a  systematic  RANS  study  of  low-Re  SD7003  airfoil  case  has 
been  performed  in  Section  4.0.  According  to  Unsteady  Reynolds  Averaged  Navier  Stokes 
(URANS)  results  presented  in  [17],  Menter’s  2-layer  Baseline  Shear  Layer  (BSL)  model  [25] 
and  a  low-Re  version  of  Jones-Launder’s  k-s  model  [12]  give  the  most  accurate  prediction  of 
Reynolds  stress.  So,  we  have  implemented  Menter’s  2-layer  BSL  model  and  Jones-Launder’s  k-e 
model  [26]  in  NASA  Langley’s  CFL3D  code  and  compare  with  the  existing  S-A  model  [27]  and 
Abid’s  low-Re  k- e  model  [13]  in  the  code.  A  simple  procedure  is  developed  to  locate  the  position 
of  separation  induced  transition  for  RANS  simulation  of  low-Re  airfoil  aerodynamics.  The 
robustness  of  the  approach  is  further  validated  with  several  other  low-Re  airfoil  cases. 

1.4  Quasi-3D  RANS  Simulation  of  Low-Re  SD7003  Airfoil  Case 

Quasi-3D  RANS  simulation  of  the  same  low-Re  airfoil  case  is  performed  with  three  planes  and 
periodic  boundary  condition  in  the  spanwise  direction.  The  objective  is  to  investigate  the  vortex 
stretching  effect  on  the  predicted  laminar  separation  bubble.  The  results  will  be  presented  in 
Section  5.0. 
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1.5  RANS  Simulation  of  University  of  Florida  5in.  Micro  Air  Vehicle  (MAV)  Wing 

The  improved  RANS  simulation  method  developed  above  will  be  applied  to  University  of 
Florida  5in.  MAV  wing.  The  computational  grids  provided  by  Prof.  Wei  Shyy  are  refined  in  the 
normal  direction.  The  results  will  be  presented  in  Section  6.0  together  with  the  laminar  solutions. 

1.6  Proper  Orthogonal  Decomposition  (POD)  Reconstruction  of  MAV  Wing  Solution 

After  obtaining  CFD  solutions  at  the  selected  training  points,  the  so-called  Proper  Orthogonal 
Decomposition  (POD)  technique  is  further  used  to  explore  the  reduced-order  modeling.  The 
results  will  be  presented  in  Section  7.0. 
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2.0  PROBLEM  STATEMENT 


Recent  advances  in  mechanical  and  electrical  system  miniaturization  have  spawned  a  growing 
interest  in  development  of  MAV’s.  These  vehicles  can  be  used  to  perform  the  military  missions 
like  battlefield  management  and  damage  assessment,  biological  and  chemical  agent  detection, 
aerial  reconnaissance,  covert  imaging  as  well  as  the  civilian  missions  like  traffic  control, 
communications  and  data  relay,  and  urban  intelligence  gathering. 

A  MAV  is  an  order  of  magnitude  smaller  than  a  conventional  one.  Accordingly  the  tip  Re  ranges 
from  about  20,000  to  500,000.  Available  in  the  literatures,  there  are  only  several  experimental 
investigations  into  the  airfoil  aerodynamics  at  such  low  Re  [1-5].  Compared  with  our 
understanding  of  high  Re  aerodynamics,  our  knowledge  of  low  Re  flows  is  very  limited  although 
there  has  been  a  long  history  of  natural  flight  studies  (e.g.,  insects  and  small  birds). 


2.1  Aerodynamics  at  Low  Re 


According  to  the  limited  studies  in  [1-5],  the  airfoil  aerodynamics  at  low  Re  is  found  quite 
different  from  our  more  familiar  aerodynamics  at  high  Re.  Due  to  large  viscous  effect,  the  flow 
at  low  Re  is  usually  dominated  by  the  formation  of  the  leading  edge  separation  bubble  and  other 
separation  phenomena,  which  degrade  the  airfoil  performance.  Figure  1  presents  the  drag  polar 
of  the  Eppler  387  airfoil,  a  well-documented  low  Re  airfoil  section,  from  [6],  A  spike  in  the 
airfoil  drag  (Cd)  is  found  for  the  lift  coefficients  (Q)  between  0.5  and  1.0  because  of  the 
formation  and  evolution  of  a  leading  edge  separation  bubble  on  the  upper  surface  of  the  airfoil. 
Therefore,  the  accurate  prediction  of  airfoil  performance  at  low  Re  requires  capturing  the  leading 
edge  separation  bubble  computationally. 
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Figure  1  -  Eppler  387  Airfoil  Lift/Drag  Polar  (taken  from  [6]) 
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The  leading  edge  separation  bubble  is  known  as  the  result  of  a  laminar  boundary  layer  separation 
due  to  a  large  adverse  pressure  gradient  after  the  point  of  minimum  pressure.  As  shown  in  Figure 
2,  after  the  separation,  the  amplification  of  velocity  disturbances  causes  the  separated  shear  layer 
transition  to  turbulence.  The  resulting  turbulent  shear  layer  extracts  more  kinetic  energy  from  the 
external  flows  and  causes  the  flow  reattached,  forming  a  separation  bubble. 


Figure  2  -  Leading  Edge  Separation  Bubble  (taken  from  [2]) 


The  transition  of  the  separated  shear  layer  from  laminar  to  turbulent  flow  is  the  key  in  the 
formation  and  evolution  of  the  leading  edge  separation  bubble.  It  determines  the  size  and  shape 
of  the  bubble  as  well  as  the  following  turbulent  boundary  layer  development.  Many  factors  affect 
the  transition  position  and  width: 

Reynolds  number 

Angle  of  attack 

Free  stream  disturbance 

Thickness  of  the  boundary  layer  at  separation 

Surface  roughness. 

The  accurate  prediction  of  this  transition  of  the  separated  shear  layer  is  really  a  challenging 
problem. 

2.2  Prediction  of  Leading  Edge  Separation  Bubble 

From  the  above  discussion,  it  is  clear  that  the  accurate  prediction  of  airfoil  performance  at  low 
Re  requires  capturing  leading  edge  separation  bubble  numerically.  As  shown  in  Figure  3,  which 
is  taken  from  [7],  without  capturing  the  leading  edge  separation  bubble,  even  the  advanced 
Navier-Stokes  approach  is  still  unable  to  accurately  predict  the  airfoil  performance  at  low  Re. 
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Figure  3  -  Results  Showing  Importance  of  Predicting  Leading  Edge  Separation  Bubble  at  Low 

Re(taken  from  [7]) 

This  situation  prompts  ZONA  Technology  to  develop  a  computational  algorithm  suitable  for 
prediction  of  the  leading  edge  separation  bubble.  Since  the  transition  of  the  separated  shear  layer 
from  laminar  to  turbulent  flow  is  the  key  in  the  formation  and  evolution  of  the  leading  edge 
separation  bubble,  the  focus  of  our  work  is  to  develop  the  capability  of  predicting  the  transition 
from  laminar  to  turbulent  flow. 
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3.0  A  FIRST  TRY  USING  LOW-Re  VERSION  OF  k-co  MODEL 


Various  computational  approaches  have  been  explored  for  prediction  of  low-Re  aerodynamics, 
ranging  from  viscous/inviscid  interactive  methods  (e.g.,  [8-11])  through  RANS  (e.g.,  [12-16]) 
and  URANS  (e.g.,  [17])  to  large-eddy  simulation  (LES)/direct  numerical  simulation  (DNS)  (e.g., 
[17-18])  methods.  Except  RANS  using  low-Re  version  of  turbulence  models  (e.g.,  [12-15])  and 
LES/DNS  (e.g.,  [17-18])  in  which  the  separation  induced  transition  occurs  naturally,  all  other 
approaches  require  an  external  mechanism  to  find  the  transition  position.  The  method  used  in 
practice  is  still  the  semi-empirical  eN  method  [19-20]  and  its  various  approximations.  In  [8],  for 
instance,  Shum  and  Marsden  used  Van  Ingen’s  shortcut  eN  method  for  fast  determination  of  the 
transition  location.  Drela  and  Giles  came  out  a  more  sophisticated  approximation,  the  so-called 
envelop  method  in  [9]  whereas  Dini  et  al.  adopted  the  table  lookup  approach  in  [10],  Stock  and 
Haase,  Yuan  et  al.  further  combined  RANS  and  URANS  with  the  eN  method  in  [16]  and  [17] 
respectively.  The  importance  of  sufficient  grid  resolution  inside  boundary  layer  has  been 
emphasized  in  [16]  for  the  accuracy  of  stability  analysis  based  on  the  Navier- Stokes  laminar 
solutions.  On  the  other  hand,  although  no  need  for  an  external  transition  mechanism,  LES/DNS 
approach  is  computationally  too  intensive  for  engineering  applications  and  the  robustness  of 
using  a  low-Re  version  of  turbulence  model  for  transition  prediction  is  questionable.  Whereas 
Wilcox  presented  some  promising  results  for  flat  plate  case  in  [14],  the  same  low-Re  version  of 
Wilcox’s  k-co  model  is  found  to  predict  the  transition  location  too  earlier  in  [15],  Furthermore, 
the  transitional  result  given  by  a  low-Re  version  of  turbulence  model  is  also  found  initial 
condition  dependent  [13], 

In  this  section,  as  a  warm-up  exercise,  we  will  first  test  the  low-Re  version  of  Wilcox’s  k-co 
model  [14]  for  the  low-Re  SD7003  airfoil  case  presented  in  [21],  Then  in  Section  4.0,  we  will 
further  develop  a  simple  method  to  locate  the  position  of  separation  induced  transition  for  RANS 
simulation  of  low-Re  airfoil  aerodynamics. 


Transition  Specific  Wilcox 's  k- co  TurbulenceModd 

According  to  [14],  the  low  Re  version  of  Wilcox’s  A:-®  turbulence  model  can  be  written  as 


with 


d  ,  d  dui  *  d  *  dk 

-(pk)  +  —(pUjk)  =  Tij- - p  pcok  +  —  [( p  +  o  nT)—] 

ot  dXj  dXj  dXj  dXj 

d  d  co  dip  2  5  ,  5©  . 

—( P«>)  +  —(pUj(o)  =  a  —  rij  — - ppco  +—[(p  +  opT  )—] 

ot  dXj  k  dXj  dXj  dXj 


(1) 


[lt  -  a  pk  /  co 


(2) 

(3) 
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(4) 


J  2  dx  : 


8u • 

+—) 

dX: 


where  t  is  time,  Xj  is  spatial  coordinate,  p  is  density,  Uj  is  velocity  component,  k  is  the  turbulent 
kinetic  energy,  co  is  specific  dissipation  rate,  %•  is  Reynolds  stress  tensor,  p  is  molecular 
viscosity,  and  pr  is  eddy  viscosity. 


In  the  form,  the  only  difference  of  the  low  Re  version  of  (1-4)  from  the  original  Wilcox’s  k-co 
turbulence  model  is  that  there  is  a  coefficient  of  a  in  front  of  pk/co  in  (3)  to  introduce  some 
low  Re  effect  with  a*  defined  as 


*  ccq  -I -Rcj*/ Rk 
1  +  Rcj/  Rk 

Other  modifications  are  all  included  in  the  closure  coefficients  of  a  and  ff: 

a  =  -fl°+V4  (a)~l 
9  1  +  ReT/Ra 

9  5/lS  +  (ReT/Rfi)4 

100  1  +  (ReT/Rp)A 


(5) 


(6) 


while  the  other  three  closure  coefficients  of  /?,  a;  a  are  the  same  as  those  in  the  original  model: 
/?  =  3/40  and  a*  =<7  =  1/2.  Here  a*0=fi/3,  a0  =  1/10,  Rp=  8,  Rk=  6,  Rco=2.1,  and 

ReT=£t. 

cop 


With  (3),  the  ru-equation  in  (1)  can  be  rewritten  as 


— ( po>)  + —( pujco)  =  ap[2a*  (Stj 
at  oxj 


1  duk  2  dUf  2  5  da)  7 

T ~z  )  —— cody  ]  — - Ppm  +  —  [(p  +  apT)—J  (7) 

3  ux  ^  3  ox  j  ox  j  ox  j 


which  is  uncoupled  from  the  /c-equation  in  (1).  Therefore,  the  k-co  model  has  a  nontrivial  laminar 
solution  for  co.,  and  we  can  start  a  computation  in  a  laminar  region  with  k= 0  in  the  boundary 
layer  and  a  small  freestream  value  of  k.  Initially  the  production  terms,  the  first  two  terms  at  the 
right  hand  side  of  both  k-  and  ^equations,  are  smaller  than  the  dissipation  term,  the  last  term  at 
the  right  hand  side  of  both  equations.  Neither  k  nor  co  is  amplified  and  the  flow  remains  laminar. 
After  a  point  where  the  critical  Re  is  achieved,  the  production  terms  start  to  exceed  the 
dissipation  term  in  the  ^-equation,  and  the  turbulent  energy  is  amplified.  At  some  point  in  this 
process,  the  eddy  viscosity  grows  rapidly  and  this  corresponds  to  the  transition  point.  After  the 
production  terms  in  the  (^-equation  further  catches  the  dissipation  term,  co  is  amplified  and 
continues  growing  until  a  balance  between  production  and  dissipation  is  achieved  in  the  k- 
equation.  Then  the  transition  from  laminar  to  turbulent  flow  is  completed  [14]. 
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Results  and  Discussion 


The  low-Re  SD7003  airfoil  case  in  [21]  is  selected  for  a  systematic  RANS  study  because  a  stable 
long  separation  bubble  is  found  over  the  upper  surface  of  the  airfoil  in  the  tests.  The  AFRL  water 
tunnel  test  conditions  are  adopted  as  the  flow  conditions  for  computations.  According  to  [21],  the 
freestream  Re  based  on  the  airfoil  chord  length  is  60,000  and  the  angle  of  attack  is  4°.  The 
freestream  Mach  number  and  turbulence  level  are  chosen  as  0.0005  and  0.08  percent 
respectively.  We  also  modify  y  and  Pr  to  fit  the  water  property. 

The  SD7003  airfoil  profile  given  by  AFRL  is  shown  in  Figure  4.  There  are  280  points  on  the 
surface.  We  have  added  another  72  points  to  extend  the  grid  to  5  chords  in  the  wake  region.  A 
hyperbolic  grid  generator  is  used  to  generate  the  computational  grid  around  the  airfoil.  There  are 
91  points  in  the  normal  direction  with  the  first  grid  spacing  of  10"5  chords  to  resolve  the  viscous 
layer. 


Figure  4  —  Computational  Grid 


The  CFD  code  used  in  this  study  is  CFL3D  v6  [22-23],  which  is  developed  and  supported  by 
NASA  Langley.  This  is  a  three-dimensional  thin-layer  RANS  solver,  using  an  implicit, 
approximately  factored,  finite  volume,  upwind  and  multigrid  algorithm.  It  employs  formally 
third-order  upwind-biased  spatial  differencing  for  the  inviscid  terms  with  flux  limiting  in  the 
presence  of  shocks.  Both  flux-difference  and  flux-vector  splitting  methods  are  available  in  the 
code.  The  flux-difference  splitting  method  of  Roe  is  employed  in  the  present  computations  for 
accurate  viscous  computations.  On  the  other  hand,  the  viscous  terms  are  discretized  with  second- 
order  central  differencing.  Since  CFL3D  is  a  compressible  CFD  code,  Weiss-Smith  low-Mach- 
number  preconditioning  technique  [24]  is  used  in  the  code  to  improve  the  convergence 
performance  for  incompressible  flow  simulation. 

Figure  5  presents  the  predicted  velocity  contours  and  streamlines  with  comparison  to  the  AFRL 
test  data.  It  is  found  that  our  laminar  result  is  very  close  to  the  test  data  except  that  the  laminar 
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result  gives  two  shorter  separation  bubbles  instead  of  a  single  long  separation  bubble  found  in 
test.  On  the  other  hand,  the  use  of  the  original  Wilcox’s  k-co  turbulence  model  gives  the  flow 
without  separation  bubble.  However,  the  use  of  the  low-Re  version  of  Wilcox’s  k-a>  turbulence 
model  improves  the  prediction  of  the  flow  pattern  given  by  the  laminar  approach.  Same  as  the 
test  data,  only  one  separation  bubble  is  predicted. 


CFD  data  (Laminar  flow ) 


x/c 


(b) 

CFD  data  (  k  -to  turbulent  model) 


Figure  5  —  Velocity  Contours  and  Streamlines  around  SD7003  Airfoil 
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4.0  RANS  SIMULATION  OF  LOW-RE  AIRFOIL  AERODYNAMICS 


Although  the  above  results  using  the  low-Re  version  of  Wilcox’s  A:-®  model  are  promising,  as 
mentioned  earlier,  the  robustness  of  using  a  low-Re  version  of  turbulence  model  for  transition 
prediction  is  questionable.  In  this  section,  we  will  develop  a  simple  method  for  RANS  simulation 
to  determine  the  position  of  separation  induced  transition  from  the  obtained  baseline  laminar 
solution,  similar  to  LES/DNS  approaches. 

Again,  the  low-Re  SD7003  airfoil  case  in  [21]  is  selected  for  a  systematic  RANS  study.  Instead 
of  using  AFRL’s  water  tunnel  conditions,  however,  the  Technische  Universitat  Braunschweig 
(TU-BS)  wind  tunnel  test  conditions  are  adopted  here  as  the  flow  conditions  for  computations. 
According  to  the  standard  atmosphere  property,  we  estimate  the  Mach  number  to  be  0.013.  The 
computational  grid  used  above  is  refined  in  the  normal  direction  with  151  points  instead  of  91 
points  and  the  first  grid  spacing  from  the  surface  is  10"6  of  the  chord  length.  The  outer  boundaries 
of  the  grid  are  further  extended  to  20  chords  in  all  directions. 

4.1  Laminar  Solution  of  Low-Re  SD7003  Airfoil  Case 

A  laminar  computation  is  performed  first.  The  predicted  surface  Cp  distributions  over  the  upper 
surface  of  the  airfoil  after  5,000,  10,000,  15,000,  20,000,  25,000,  30,000,  and  35,000  iterations 
are  presented  in  Figure  6.  It  is  found  that  after  10,000  iterations,  the  obtained  surface  pressure 
coefficient  (Cp)  distribution  is  still  like  the  one  without  the  bubble  perturbation.  By  that  time,  as 
shown  in  the  convergence  history  presented  in  Figure  7,  the  numerical  residual  has  already 
dropped  by  three  orders  of  magnitude.  So,  most  people  may  stop  the  computation  there.  They 
thought  that  a  fully  converged  solution  is  obtained  and  without  any  external  transition 
mechanism,  the  existing  CFD  approach  is  not  able  to  capture  the  laminar  separation  bubble. 
However,  if  one  is  not  fooled  by  the  misleading  convergence  history  and  keeps  running  the 
computation,  after  15,000  iterations,  as  shown  in  Figure  6,  a  surface  Cp  distribution  with  the 
bubble  perturbation  starts  to  be  captured  in  the  computations,  with  a  Cp  plateau  between  the 
separation  and  transition.  After  25,000  iterations,  the  predicted  surface  Cp  distribution  becomes 
quite  stable.  This  is  different  from  Pauley  et  al.’s  finding  in  [28],  where  a  two-dimensional 
laminar  separation  solution  is  found  quite  unstable. 
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x/c 


Figure  6  -  Predicted  Surface  Cp  Distributions  on  the  Upper  Surface  of  SD7003  Airfoil 


Also  shown  in  Figure  6  is  the  3D  LES  result  from  [17].  It  is  found  that  our  converged  laminar 
solution  agrees  with  the  LES  solution  relatively  well  up  to  the  transition  point.  After  the 
transition  point,  there  is  a  much  larger  discrepancy  between  the  two  solutions.  Most  notably, 
different  from  the  LES  result,  our  converged  laminar  solution  has  several  artificial  humps  in  the 
surface  Cp  distribution  after  the  transition  point.  Looking  from  DNS  point  of  view,  this  indicates 
that  our  employed  numerical  method  and  mesh  are  able  to  resolve  the  simulated  laminar  flow 
quite  well  but  still  unable  to  resolve  the  turbulent  flow  yet.  A  proper  turbulence  model  is 
required  for  simulation  of  the  flow  after  the  transition  point. 


Figure  7  -  Convergence  History 
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It  is  noteworthy  that  after  9,000  iterations,  as  shown  in  Figure  7,  the  numerical  residual  quickly 
jumps  by  two  orders  of  magnitude  and  oscillates  around  that  level.  As  will  be  shown  later,  this  is 
solely  because  we  do  not  turn  on  a  turbulence  model  from  the  transition  point  in  the 
computation.  Based  on  the  surface  Cp  distribution  shown  in  Figure  6  and  the  velocity  fields 
shown  in  the  following  Figure  8,  the  big  jump  of  numerical  residual  after  9,000  iterations 
indicates  the  transition  of  the  numerical  solution  from  one  without  the  bubble  perturbation  to  the 
one  with  the  bubble  perturbation. 

Figure  8  further  presents  the  predicted  velocity  fields,  including  some  streamlines,  after  5,000, 
10,000,  15,000,  20,000,  25,000,  30,000,  and  35,000  iterations  and  compares  with  AFRL  test  data 
in  [21].  It  is  found  that  the  flowfield  solution  obtained  after  5,000  iterations  has  no  separation  at 
all.  By  10,000  iterations,  the  trailing-edge  separation  is  captured  and  by  15,000  iterations,  the 
laminar  separation  bubble  further  starts  to  be  captured.  After  25,000  iterations,  the  predicted  flow 
patterns  are  quite  stable.  All  these  observations  are  consistent  with  the  surface  Cp  surface 
distributions  shown  in  Figure  6.  Compared  with  AFRL  test  data,  the  converged  laminar  solution 
of  the  separation  bubble  is  much  longer  and  thicker  because  there  are  three  vortical  structures 
inside  the  separation  bubble  instead  of  one  found  in  the  test.  This  also  explains  those  artificial 
humps  found  in  the  converged  surface  Cp  distribution  shown  in  Figure  6.  Apparently,  the  latter 
two  vortical  structures  are  due  to  the  lack  of  transition  to  turbulence  in  the  computation  and  they 
are  not  physical. 
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(a)  after  5,000  iterations  (b)  after  10,000  iterations 


(c)  after  15,000  iterations 


(d)  after  20,000  iterations 


0  0.5 

x 

(g)  after  35,000  iterations 


(f)  after  30,000  iterations 


-0.1 

-0.2 

-0.3 


(h)  AFRL  test  data 


Figure  8  -  Velocity  Field  and  Streamlines 


13 


Besides  insufficient  iterations,  another  possible  reason  for  people  to  miss  capturing  a  stable 
laminar  separation  bubble  in  computation  is  the  insufficient  grid  resolution.  For  a  typical  viscous 
computation  like  the  current  one,  a  mesh  with  91  points  in  the  normal  direction  and  the  first  grid 
spacing  of  10~5  of  the  chord  length  from  the  surface  is  usually  sufficient  to  resolve  the  boundary 
layer  flows.  However,  as  shown  in  Figure  9,  the  surface  Cp  distribution  predicted  on  this  coarser 
mesh  is  very  unstable,  the  same  as  Pauley  et  al.’s  finding  in  [28].  Apparently  although  the 
resolution  of  this  grid  is  sufficient  to  resolve  the  boundary  layer,  it  is  still  not  sufficient  to  resolve 
the  laminar  flow  separation. 


Figure  9  -  Predicted  Surface  Cp  Distribution  on  a  Coarser  Mesh  with  91  Points  in  Normal 

Direction 


4.2  Specification  of  Transition  Point 

A  key  for  prediction  of  laminar  separation  bubble  is  how  to  specify  the  separation  induced 
transition  point.  The  existing  transition  prediction  methods  range  from  simple  empirical 
correlations  through  semi-empirical  methods  like  eN  method  to  direct  numerical  simulations.  The 
eN  method,  based  on  linear  stability  analysis  and  the  empirical  input  of  N  at  transition,  is  the 
most  popular  in  practice.  In  [16],  Stock  and  Haase  investigated  the  feasibility  of  coupling  eN 
method  to  RANS  computation.  It  is  found  that  in  order  to  produce  reliable  laminar  data  for  the 
stability  analysis,  a  sufficiently  large  and  constant  number  of  grid  points  are  required  inside  the 
viscous  layer.  This  requirement,  combined  with  the  stability  analysis,  render  the  resulting  RANS 
computation  much  more  intensive.  In  the  following,  we  will  present  a  much  simpler  method  to 
locate  the  transition  point. 

Let  us  blow  up  Figure  8(g)  near  the  upper  surface  of  the  airfoil.  As  shown  in  Figure  10,  the  flow 
next  to  the  surface  reverses  the  direction  around  0.18c,  indicating  the  separation  of  laminar  flow. 
From  the  separation  point,  a  dividing  line  forms.  Above  this  line,  the  flow  moves  downstream 
and  below  this  line,  the  flow  moves  upstream,  resulting  in  a  strong  free  shear  layer  there.  This 
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shear  layer  starts  to  become  unstable  from  0.4c  and  eventually  rolls  up  into  a  vortex.  At  the  same 
time,  the  reversed  flow  next  to  the  upper  surface  reverses  the  direction  again  and  moves 
downstream  whereas  in  reality  this  should  happen  until  the  reattachment.  Both  phenomena  are 
not  physical  and  can  be  used  as  the  indications  of  the  fact  that  a  laminar  solution  is  unsustainable 
from  there  and  the  transition  to  turbulence  should  take  place. 


Figure  10  -  Velocity  Fields  Adjacent  to  the  Upper  Surface  of  SD7003  Airfoil 

To  quantify  the  location  of  the  transition  point,  the  distribution  of  the  tangential  Mach  numbers 
at  the  grid  points  next  to  the  upper  surface  of  SD7003  airfoil  is  further  presented  in  Figure  1 1.  It 
is  found  that  the  laminar  solution  is  very  smooth  before  the  second  flow  reverse  but  becomes 
quite  unstable  after  the  second  flow  reverse.  This  further  proves  that  the  transition  to  turbulence 
should  take  place  where  the  laminar  solution  next  to  the  upper  surface  reverses  the  direction  for 
the  second  time.  According  to  Figure  1 1,  the  transition  point  for  this  low-Re  SD7003  airfoil  case 
is  at  0.416c  whereas  the  laminar  separation  point  is  at  0.176c. 


Figure  11  -  Distribution  of  Tangential  Mach  Numbers  next  to  the  Upper  Surface  of  SD7003  Airfoil 

Also  included  in  Figure  11  are  the  “laminar”  solutions  from  RANS  simulations  with  zero 
production  terms  in  the  selected  turbulence  models.  The  use  of  such  RANS  solutions  as  the 
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“laminar”  solutions  has  two  advantages.  First,  this  will  introduce  less  perturbation  to  the 
subsequent  RANS  computations  in  which  the  production  terms  in  the  turbulence  models  are 
turned  on  after  the  transition  point.  Second  but  more  importantly,  the  effect  of  the  freestream 
turbulence  level  on  the  transition  point  can  be  taken  into  account  through  the  boundary 
conditions.  It  is  found  that  the  two  solutions  given  by  the  laminar  computation  and  the  RANS 
simulation  with  Spalart-Allmaras  (S-A)  one-equation  model  [27]  are  almost  indistinguishable  up 
to  the  transition  point.  On  the  other  hand,  the  solution  given  by  the  RANS  simulation  with  Jones- 
Launder’s  (J-L)  k-e  model  [26]  deviates  from  the  laminar  solution  the  most.  The  transition  point 
given  by  the  RANS  simulation  with  Menter’s  2-layer  BSL  model  [25]  is  one  grid  point  later  than 
the  laminar  solution  and  the  S-A  solution  whereas  J-L  model  delays  the  predicted  transition 
position  by  another  grid  point  from  the  BSL  solution.  All  these  findings  are  further  supported  by 
the  surface  Cp  distributions  shown  in  Figure  12.  Among  all  RANS  solutions,  the  S-A  solutions 
are  the  closet  to  the  laminar  solutions  in  Figure  6.  On  the  other  hand,  the  initial  discrepancy 
between  the  RANS  solutions  given  by  the  two-equation  models  and  the  laminar  solutions  is  quite 
large.  However,  the  differences  between  the  solutions  up  to  the  transition  point  become  smaller 
with  iterations.  From  the  transition  point,  moreover,  the  solutions  given  by  RANS  with  the  two- 
equation  models  are  found  less  stable  than  the  laminar  solutions. 
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(a)  S-A  model 


0  0.1  0.2  0.3  0.4  0.5  0.6  0.7  0.8  0.9  1 

x/c 

(b)  Menter  2-layer  BSL  model 


x/c 


(c)  Jones-Launder  k-s  model 


Figure  12  -  Surface  Cp  Distributions  Predicted  by  RANS  with  Zero  Production  Terms 
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4.3  RANS  Solution  of  Low-Re  SD7003  Airfoil  Case 


After  locating  the  transition  point,  RANS  simulations  are  performed  with  zero  production  terms 
in  the  selected  turbulence  model  before  the  transition  point  and  nonzero  production  terms  after 
the  transition  point.  To  compare  the  convergence  performance  with  the  above  laminar 
computation,  we  start  RANS  computation  from  the  uniform  freestream  conditions  again  rather 
than  restart  the  computation  from  the  above  RANS  solutions  with  zero  production  terms  in  the 
selected  turbulence  model.  It  is  found  in  Figure  13  that  initially  the  numerical  residual  of  RANS 
simulation  with  S-A  model  is  almost  the  same  as  the  laminar  solution.  However,  different  from 
the  laminar  solution,  the  numerical  residual  of  RANS  simulation  does  not  jump  back  after  9,000 
iterations.  Instead  the  numerical  residual  drops  monotonically.  This  is  due  to  turn-on  of  the 
production  term  in  the  selected  S-A  model  from  the  transition  point. 


Figure  13  -  Convergence  History 

The  predicted  surface  Cp  distributions  over  the  upper  surface  of  the  airfoil  by  RANS  simulations 
with  different  turbulence  models  are  presented  in  Figure  14  along  with  the  3D  LES  result  in  [17]. 
It  is  found  that  the  converged  S-A  solution  is  almost  indistinguishable  from  the  3D  LES  result. 
On  the  other  hand,  the  solutions  given  by  the  two-equation  models  are  quite  different  from  the 
3D  LES  result.  It  is  interesting  to  note  that  different  from  Horton’s  bubble  model,  the  S-A 
solution  does  not  immediately  accelerate  the  deceleration  process  after  the  transition  point. 
Instead  the  predicted  flow  keeps  the  same  deceleration  rate  for  a  while  and  starts  to  decelerate 
quickly  until  after  0.532c.  On  the  other  hand,  consistent  with  Horton’s  model,  the  converged 
BSL  solution  speeds  up  the  deceleration  immediately  after  the  transition  point.  The  predicted 
deceleration  rate  is  also  larger  than  the  3D  LES  result,  indicating  the  possibility  of  overpredicting 
Reynolds  stress  there.  Furthermore,  the  converged  J-L  solution  is  found  between  the  above  two 
RANS  solutions.  The  deceleration  after  the  transition  point  speeds  up  before  the  S-A  solution  but 
after  the  BSL  solution.  The  predicted  deceleration  rate  is  smaller  than  those  given  by  the  other 
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approaches,  indicating  the  possibility  of  underpredicting  Reynolds  stress  there.  Also  included  in 
Figure  14  is  the  result  given  by  RANS  simulation  with  Abid’s  low-Re  version  of  k-e  model 
(Abid)  [13]  without  specification  of  the  transition  point.  It  is  found  that  this  low-Re  version  of  k- 
s  model  does  capture  the  laminar  separation  bubble  automatically  but  is  unable  to  completely 
remove  those  humps  found  in  the  above  laminar  Cp  distribution  associated  with  the  artificial 
vortical  structures  inside  the  boundary  layer.  Therefore,  the  flow  reattaches  the  surface  much 
later  than  the  3D  LES  result  and  the  separation  bubble  is  unreasonably  long  and  thick. 


Figure  15  presents  the  predicted  Reynolds  stress  distributions  by  RANS  with  the  three  turbulence 
models  in  comparison  with  AFRL  test  data  in  [21].  It  is  surprising  to  see  that  all  RANS 
computations,  especially  Menter’s  2-layer  BSL  model,  underpredict  the  maximum  absolute  value 
of  Reynolds  stress.  This  is  inconsistent  with  the  predicted  surface  Cp  distributions  shown  in 
Figure  14.  A  possible  reason  for  that  is  we  compare  our  RANS  results  with  AFRL  water  tunnel 
test  data  instead  of  TU-BS  wind  tunnel  test  data.  According  to  [17],  the  maximum  absolute  value 
of  Reynolds  stress  in  AFRL  water  tunnel  test  data  is  slightly  larger  than  TU-BS  wind  tunnel  test 
data.  On  the  other  hand,  the  starting  position  and  shape  of  the  predicted  Reynolds  stress  pocket 
by  S-A  and  J-L  models  agree  with  the  test  data  well  whereas  the  Reynolds  stress  pocket  given  by 
BSL  model  starts  too  early.  Especially  the  predicted  Reynolds  stress  by  BSL  model  jumps  to  its 
maximum  absolute  value  too  quickly  after  the  transition  point.  This  is  why  the  BSL  surface  Cp 
solution  has  a  turning  point  immediately  after  the  transition.  Consistent  with  the  surface  Cp 
distributions  shown  in  Figure  14,  BSL  model  produces  the  largest  maximum  absolute  value  of 
Reynolds  stress  among  all  RANS  solutions  whereas  J-L  model  gives  the  smallest  one. 
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Figure  15  -  Reynolds-stress  Distributions 

To  further  compare  the  predicted  laminar  separation  bubble  with  AFRL  test  data,  Figure  16 
presents  the  velocity  contours  and  streamlines  around  SD7003  airfoil.  It  is  found  that  the 
separation  bubble  predicted  by  S-A  model  is  the  closest  to  the  test  data,  including  both  the  length 
and  thickness,  whereas  those  given  by  the  two-equation  models  are  too  short  and  too  thin. 
Especially  the  one  from  BSL  model  is  almost  negligible. 
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Figure  16  -  Velocity  Contours  and  Streamlines 
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These  facts  can  be  more  clearly  seen  in  the  distributions  of  tangential  Mach  numbers  next  to  the 
upper  surface  of  airfoil  shown  in  Figure  17.  It  is  found  that  after  turning  on  the  production  terms 
in  the  selected  turbulence  model  after  the  transition  point,  the  laminar  separation  point  can  be 
modified  to  a  significantly  later  position.  Especially  for  the  BSL  solution,  there  are  only  two  grid 
points  in  the  reversed  flow  state  due  to  overprediction  of  Reynolds  stress. 


Figure  17  -  Distribution  of  Tangential  Mach  Numbers  next  to  the  Upper  Surface  of  SD7003  Airfoil 

Table  1  further  presents  a  quantitative  comparison  between  the  numerical  separation  bubble  and 
the  existing  test  data  in  [21].  Compared  with  the  test  data  and  3D  LES  result,  S-A  model  is  found 
to  predict  a  longer  bubble  whereas  the  two-equation  models  give  a  shorter  bubble.  Again  the 
separation  bubble  predicted  by  2-layer  BSL  model  is  almost  negligible. 


Table  1  -  Comparison  of  Numerical  and  Experimental  Results 


xs/c 

xt/c 

xr/c 

Test  (TUBS) 

0.30 

0.55 

0.62 

Test  (AFRL) 

0.18 

0.47 

0.58 

3D  LES 

0.25 

0.49 

0.60 

S-A 

0.22 

0.42 

0.67 

BSL 

0.40 

0.42 

0.42 

J-L 

0.30 

0.43 

0.51 

4.4  Application  to  Other  Low-Re  Airfoil  Cases 

To  examine  the  robustness  of  the  developed  method  for  prediction  of  low-Re  airfoil 
aerodynamics,  we  will  further  investigate  several  other  low-Re  airfoil  cases.  This  time  RANS 
computations  with  zero  production  term  in  the  selected  S-A  model  are  performed  first  to  provide 
the  baseline  “laminar”  solution.  After  determination  of  the  transition  point  using  the  method 
developed  above,  RANS  computations  with  the  complete  S-A  model  after  the  transition  point  are 
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restarted  from  the  above  RANS  “laminar”  solutions  instead  of  the  uniform  freestream  conditions 
for  less  numerical  iterations  before  a  fully  converged  RANS  solution  is  achieved. 

Epplet'387  Airfoil 

The  first  low-Re  airfoil  case  investigated  is  the  Eppler387  airfoil  case  in  [3].  It  is  noteworthy  that 
this  test  model  has  a  blunted  trailing  edge  with  the  thickness  of  0.167  percent  of  the  chord  length. 
Our  numerical  results  indicate  that  whether  or  not  modeling  this  blunt  trailing  edge  has  a  big 
impact  on  the  predicted  Cp  suction  peak  value,  especially  for  higher  angle  of  attack.  Figure  18 
presents  the  two  types  of  meshes  we  have  used  in  computation.  The  first  one  is  a  C-mesh  without 
modeling  the  blunt  trailing  edge.  There  are  301  points  in  the  wraparound  direction  with  201 
points  on  the  body,  and  151  points  in  the  normal  direction  with  the  first  grid  spacing  from  the 
surface  as  10"6  of  the  chord  length.  The  second  one  is  an  O-mesh  with  modeling  the  blunt  trailing 
edge.  There  are  25 1  points  in  the  wraparound  direction  with  1 1  points  to  cover  the  blunt  trailing 
edge,  and  111  points  in  the  normal  direction  with  the  first  grid  spacing  from  the  surface  as  10"6  of 
the  chord  length. 


Figure  18  -  Computational  Grids  Used  for  Eppler387  Airfoil  Cases 

There  are  many  test  conditions  in  [3].  The  first  one  investigated  here  is  the  case  with  the 
freestream  Mach  number  of  0.08,  and  the  Re  of  100,000,  and  the  angle  of  attack  of  7°.  Figure  19 
presents  the  predicted  surface  Cp  distribution  on  the  two  computational  grids  in  comparison  with 
the  test  data  from  [3].  It  is  found  that  without  modeling  the  blunt  trailing  edge,  the  predicted  peak 
value  of  the  surface  Cp  suction  on  the  C-mesh  is  much  lower  than  the  test  data.  On  the  other 
hand,  with  modeling  the  blunt  trailing  edge,  the  predicted  surface  Cp  distribution  on  the  O-mesh 
is  nearly  the  same  as  the  test  data. 
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Figure  19-Surface  Cp  Distributions  for  Eppler387  Airfoil  Case 


Next,  we  will  further  investigate  the  most  difficult  cases  in  [3]  with  the  Re  of  60,000.  At  this  Re, 
the  flow  is  found  not  stable  in  the  tests,  making  an  accurate  prediction  very  difficult.  In  Figure 
20,  we  present  the  laminar  solutions  of  Cl  versus  a  on  both  C-mesh  and  O-mesh  and  the  RANS 
solution  using  the  S-A  model  and  transition  mechanism  on  C-mesh  in  comparison  with  the  test 
data  in  [3].  As  expected,  the  predicted  Cl  values  on  C-mesh  with  sharp  trailing  edge  are  smaller 
than  those  predicted  on  O-mesh  with  blunted  trailing  edge.  This  is  because  as  shown  in  Figure 
19,  the  former  underpredicts  the  Cp  suction  peak  value. 


Figure  20-Lift  Curve  for  Eppler387  Airfoil  Case 
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LA203A  Airfoil 


The  second  investigated  case  is  the  low-Re  LA203A  airfoil  case  in  [29].  Again  an  O-type  mesh 
is  used  for  this  case.  As  shown  in  Figure  21,  there  are  161  points  in  the  wraparound  direction  and 
121  points  in  the  normal  direction  with  the  first  grid  spacing  from  the  surface  as  10"6  of  the  chord 
length. 


Figure  21-Computational  Grid  Used  for  LA203A  Airfoil  Case 


The  case  investigated  has  the  freestream  Mach  number  of  0.1,  and  the  Re  of  250,000.  The  angle 
of  attack  is  4°.  Figure  22  presents  the  predicted  laminar  and  turbulent  solutions  of  the  surface  Cp 
distribution.  It  is  found  that  the  obtained  baseline  laminar  solution  is  already  very  close  to  the 
test  data  except  that  the  Cp  plateau  is  shorter  compared  with  the  test  data.  After  turning  on  the 
production  term  in  the  S-A  model  after  the  transition  point,  the  obtained  solution  is  further 
improved. 
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Figure  22-Surface  Cp  Distributions  for  La203a  Airfoil  Case 


LNV109A  Airfoil 

The  third  case  considered  is  the  low-Re  LNV109A  airfoil  case  in  [29],  A  C-type  mesh  is  used  for 
this  case.  As  shown  in  Figure  23,  there  are  281  points  in  the  wraparound  direction  with  181 
points  on  the  body  and  121  points  in  the  normal  direction  with  the  first  grid  spacing  from  the 
surface  as  10'6  of  the  chord  length.  The  outer  boundary  of  the  grid  is  extended  to  20  chords  in  all 
directions. 


Figure  23-Computational  Grid  Used  for  LNV109A  Airfoil  Case 
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Figure  24  presents  the  predicted  laminar  and  turbulent  solutions  of  the  surface  Cp  distribution  in 
comparison  with  the  test  data  in  [29]  for  the  case  with  the  freestream  Mach  number  of  0.1,  and 
the  Re  of  375,000,  and  the  angle  of  attack  of  4°.  It  is  found  that  after  turning  on  the  production 
term  in  the  S-A  model  after  the  transition  point,  those  numerical  humps  in  the  laminar  solution  of 
the  surface  Cp  distribution  after  the  transition  point  are  removed.  Otherwise  the  laminar  and 
RANS  solutions  are  very  similar.  Both  have  shorter  Cp  plateau  than  the  test  data.  The 
importance  of  the  baseline  laminar  solution  on  the  accuracy  of  the  final  solution  is  clearly  seen. 


Figure  24-Surface  Cp  Distributions  for  LNV109A  Airfoil  Case 


FX63-1 37  Airfoil 

The  next  case  considered  is  the  low-Re  FX63-137  airfoil  case  in  [30].  A  C-type  mesh  is  used  for 
this  case.  As  shown  in  Figure  25,  there  are  291  points  in  the  wraparound  direction  with  191 
points  on  the  body  and  121  points  in  the  normal  direction  with  the  first  grid  spacing  from  the 
surface  as  10"6  of  the  chord  length.  The  outer  boundary  of  the  grid  is  extended  to  20  chords  in  all 
directions. 
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Figure  25-Computational  Grid  Used  for  FX63-137  Airfoil  Case 


Figure  26  presents  the  predicted  laminar  and  turbulent  solutions  of  the  surface  Cp  distribution  in 
comparison  with  the  test  data  in  [30]  for  the  case  with  the  freestream  Mach  number  of  0.022,  and 
the  Re  of  100,000,  and  the  angle  of  attack  of  7°.  It  is  found  that  there  are  some  numerical 
oscillations  after  the  suction  peak  in  the  final  surface  Cp  distribution,  which  also  exist  in  the 
laminar  solution.  Apparently  this  problem  is  related  to  either  the  smoothness  of  the  surface 
geometry  definition  we  have  got  or  the  robustness  of  the  applied  low  Mach  number 
preconditioning  method  but  has  nothing  to  do  with  the  developed  transition  mechanism  in  this 
Phase  I  work.  After  turning  on  the  production  term  in  the  S-A  model  after  the  transition  point, 
the  obtained  numerical  solution  is  improved  but  still  slightly  different  from  the  test  data  in  [30], 
A  possible  reason  for  this  is  that  we  have  not  digitized  the  test  data  accurately. 
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Figure  26-Surface  Cp  Distributions  for  FX63-137  Airfoil  Case 
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5.0  QUASI-3D  RANS  SIMULATION  OF  LOW-RE  SD7003  AIRFOIL  CASE 


Next,  we  will  investigate  the  above  low-Re  SD7003  airfoil  case  again  using  a  3-D  grid  with  3 
planes  in  the  spanwise  direction  and  using  a  periodic  boundary  condition  along  the  spanwise 
direction.  Each  C-mesh  plane  is  the  same  as  the  2-D  mesh  used  before  with  386  points  along  the 
wraparound  direction,  and  280  points  on  the  airfoil  surface,  and  151  points  in  the  normal 
direction  with  the  first  grid  spacing  from  the  surface  as  10~6  of  the  chord  length.  The  grid  spacing 
in  the  spanwise  direction  is  0.4  chords.  Our  objective  is  to  investigate  the  effect  of  vortex 
stretching  on  the  laminar  separation  bubble. 


Figure  27  -  Computational  Grid 

A  laminar  computation  is  performed  first.  The  predicted  surface  Cp  distribution  over  the  upper 
surface  of  the  airfoil  at  each  spanwise  station  after  35,000  iterations  is  presented  in  Figure  28 
together  with  the  2-D  RANS  solution  and  3-D  LES  result.  It  is  found  that  there  is  no  visible 
difference  between  the  obtained  surface  Cp  distributions  at  each  spanwise  station,  indicating  no 
3-D  effect  captured  in  the  RANS  computation.  On  the  other  hand,  we  do  see  that  the  quasi-3D 
solution  is  different  from  the  2-D  solution,  especially  those  artificial  humps  in  the  predicted 
surface  Cp  distributions.  As  mentioned  in  the  last  section,  the  artificial  humps  in  the  surface  Cp 
distributions  result  from  the  artificial  vortical  structures.  The  difference  between  the  quasi-3D 
solution  and  the  2-D  solution  indicates  that  the  predicted  vortical  structures  in  the  quasi-3D 
computation  are  different  from  those  in  the  2-D  computation  due  to  vortex  stretching  effect. 
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Figure  28  -  Laminar  Surface  Cp  Distributions  on  the  Upper  Surface  of  SD7003  Airfoil 

Next,  RANS  computation  is  performed  with  the  numerical  method  developed  in  the  last  section 
using  the  S-A  turbulence  model.  The  resulting  surface  Cp  distributions  over  the  upper  surface  of 
the  airfoil  are  presented  in  Figure  29  together  with  the  2-D  RANS  result  and  3-D  LES  result. 
This  time  we  find  that  the  quasi-3D  solution  is  almost  the  same  as  the  2-D  solution  even  though 
the  3-D  transition  location,  0.407c,  is  slightly  earlier  than  the  2-D  transition  location,  0.416c. 
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Figure  29  -  S-A  RANS  Surface  Cp  Distributions  on  the  Upper  Surface  of  SD7003  Airfoil 


31 


6.0  RANS  SIMULATION  OF  UNIVERSITY  OF  FLORIDA  5  IN.  MAY  WING 


The  real  MAV  wing  configuration  we  select  for  investigation  in  Phase  I  is  University  of  Florida 
(UF)  5in.  MAV  wing,  as  shown  in  Figure  30.  This  membrane  wing  has  the  root  chord  length  of 
4.09  in.  and  the  wing  span  of  5  in.  The  wing  area  is  17.1 1  in2.  In  Phase  I,  we  only  consider  the 
rigid  wing  case. 


6”  Maximum  Dimension 
Carbon  Fiber  Skeleton 
Latex  Rubber  Skin 
52  grams,  1.8  oz. 

Hand  Launched 


Speed  12-30  mph 


Color  video  camera/transmitter 
^  Maxon  electric  motor 
Lithium  polymer  batteries 
Miniature  RC  receiver 
2  miniature  rotary  servos 
Fully  proportional  speed  controller 


Figure  30-UF  MAV:  Flexible  Latex  Rubber  Wing 


The  computational  grid  around  this  wing  configuration  is  provided  by  Prof.  Wei  Shyy.  The  grid 
was  created  using  ANSYS  ICEM  CFD  grid  generation  software  based  on  the  wing  geometry 
exported  from  ProE.  The  multiblock  structured  grid  has  18  blocks  and  210,000  grid  points.  The 
wing  surface  is  covered  by  41  by  21  points.  The  top  and  bottom  outer  boundaries  are  at  6c 
distance  from  the  wing  whereas  the  outflow  boundary  is  at  1  lc  distance  from  the  wing.  Due  to 
flow  symmetry,  only  half  of  the  wing  is  modeled. 


Figure  31-Computational  Grid  Blocks  Layout 
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Based  on  our  experience  for  SD7003  airfoil  case,  the  resolution  of  the  above  grid  is  not  sufficient 
for  accurate  simulation  of  laminar  flow  separation.  So,  we  have  refined  the  grids  in  the  normal 
direction  by  reducing  the  first  grid  spacing  from  the  surface  to  10'6  of  root  chord  length  and 
increasing  the  number  of  grid  points  in  the  normal  direction  from  25  to  99  for  each  block. 

Again  a  laminar  computation  is  performed  first.  The  flow  conditions  are:  M,y,  =  0.0294  and  Re  = 
90,000  based  on  the  root  chord  length.  To  construct  CFD  solution  database  for  POD  analysis,  a 
series  of  computations  are  performed  at  the  selected  training  points:  a  =  -5°,  0°,  5°,  ...,  50°.  The 
predicted  laminar  solution  of  surface  density  contours  and  entropy  contours  in  5  cutting  planes 
are  presented  in  the  left-hand  side  of  Figure  32.  Also  presented  in  the  right-hand  side  of  Figure 
32  are  the  RANS  solutions  using  the  S-A  turbulence  model  and  the  above  developed  transition 
mechanism.  It  is  found  that  the  surface  density  contours  predicted  by  the  two  approaches  are 
quite  different,  especially  for  higher  angles  of  attack,  where  the  turbulent  flow  regions  are  larger. 
On  the  other  hand,  the  predicted  tip  vortex  structures  by  the  two  approaches  are  quite  similar. 
This  is  because  we  assume  that  outboard  of  the  wing  tips,  the  flow  is  always  laminar. 
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(e)  a  =  35° 


Figure  32-Flow  Patterns  around  University  of  Florida  5  in.  MAY  Wing 
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7.0  POD  RECONSTRUCTION  OF  MAY  WING  SOLUTION 


A  Navier-Stokes  simulation  is  far  too  computationally  expensive  to  be  used  for  MAV  stability 
and  control  analysis.  Therefore,  ZONA  Technology  proposes  to  use  POD  to  further  distill  the 
computed  data  from  Navier-Stokes  simulations  into  reduced  order  models  that  can  serve  as  the 
basis  for  further  MAV  stability  and  control  analysis. 

POD  is  a  technique  for  extraction  of  a  set  of  basis  functions  of  the  spatial  variables  that  have 
maximum  energy  content,  and  thereby  the  capture  of  the  principal  structures,  allowing  for 
construction  of  a  model  of  reduced  dimension  to  approximate  the  original  ensemble. 

Let  us  assume  the  solution  to  be  a  series  of  steady  state  variations  of  the  angle  of  attack  whereby 
the  angle  of  attack  is  treated  just  like  time  in  conventional  POD  approaches  for  time  dependent 
fluid  flows.  Therefore,  for  a  given  Mach  number  and  Re, 

K 

fix,  a)=  1 4  (a)yrk  (x)  (8) 

k= 1 

where  f(x,a )  is  the  flow  variable,  </>k(a)  is  the  scalar  coefficient  with  a  being  the  axis  on 
which  the  flow  field  evolves,  and  y/k  (x)  is  the  basis  function.  The  same  concept  holds  for  angle 
of  attack,  Mach  number,  and  Re  as  the  axis  to  trace  the  variations  in  MAV  flow  fields. 

The  choice  of  y/k(x)  in  Eq.  (8)  is  not  unique.  To  facilitate  reduced-order  modeling,  POD 
searches  for  the  optimal  orthogonal  basis  functions  y/k{x)  in  the  sense  that  the  linear 
combination  of  the  first  few  basis  functions  gives  an  optimal  convergence  path  for  the  series  in 
Eq.  (8),  thus  justifying  its  truncation  well  before  all  the  K  combinations.  According  to  [31],  this 
can  be  achieved  by  a  singular  value  decomposition  of  the  solution  matrix  S  with  each  column 
representing  a  set  of  CFD  solutions  for  each  angle  of  attack.  As  a  result,  the  above  POD  problem 
becomes  a  search  for  K  eigenvectors  corresponding  to  the  K  largest  eigenvalues  of  the  matrix 
SST.  The  Air  Force’s  SNAPMAN2  software  [32]  is  used  to  perform  the  POD  analysis. 

After  obtaining  those  optimal  basis  functions,  the  flow  variable  in  (8)  can  be  represented  by 

fix, «)  =  I  <4  (aWk  ix)  (9) 

k= 1 

The  rapid  convergence  of  the  norm  of  E^(x,a)  =  f(x,a)-fK  (x,a)  with  K  allows  for  a 
reduced-order  modeling,  where 

fK(x,a)=  ZAiaWkix)  (10) 

k= 1 

and  K  «K .  The  percentage  of  the  flow  energy  captured  by  a  reduced-order  model  can  be 
defined  as 

%  =  £4  (U) 

k= 1  k= 1 

where  Xk  is  the  eigenvalue  of  the  matrix  SST. 
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Results  and  Discussion 


In  this  section,  the  results  from  our  POD  analysis  of  the  high-level  CFD  solution  database  for 
University  of  Florida  5  in.  MAV  wing  are  presented.  It  is  found  that  for  higher  angles  of  attack, 
an  accurate  POD  reconstruction  requires  more  modes. 


CFD  solution 


POD  reconstruction  with  4  modds 


(a)  a  =  5° 


CFD  solution  POD  reconstruction  with  4  mod6s 


(b)  a  =  25° 


CFD  solution 


POD  reconstruction  with  4  modek 


(c)  a  =  50° 


Figure  33-POD  Reconstructions  of  Surface  Cp  Distribution 
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8.0  PHASE  II  PLAN 


8.1  Phase  II  Objectives 

The  overall  objective  is  to  establish  an  enhanced  aerodynamic  methodology  to  handle  a  low-Re 
and/or  a  gusty  flight  regime  encountered  by  a  membrane-wing  MAV  for  its  control  and 
performance.  The  specific  technical  objectives  on  Phase  II  include: 

a.  Establish  an  expedient  computational  method  in  low-Re  aerodynamics  for  the  performance 
and  control  of  a  membrane  MAV  (mMAV)  wing  in  a  gusty  and  a  gust-free  environment.  The 
methodology  should  be  able  to  predict  airload/response  due  to  gust  on  a  rigid/flexible  mMA  V 
wing  with  occurring  separation  bubble  and  flow  transition. 

b.  Experimental  verifications  of  the  computational-aerodynamic  solutions  of  an  airfoil  and  a 
rigid  mMA  V  wing  in  a  gusty  and/or  a  gust-free  environment. 

c.  Using  the  enhanced  aerodynamic  method  developed  to  obtain  time-domain  forces  and 
moments  for  the  rigid/flexible  mMAV  thus  providing  inputs  for  a  3DOF ZSimulink 
demonstration  of  the  mMAV  control/  performance. 

8.2  Phase  II  Tasks 

Six  major  tasks  will  be  performed  in  Phase  II  which  include. 

(1)  Rigid  mMAV  wing  aerodynamics,  (2)  Flexible  mMAV  wing  aero  elasticity,  (3)  Gust 
Considerations,  (4)  Perform  Proper  Orthogonal  Decomposition/Response  Surface  Methodology 
(POD/RSM)  of  the  Unsteady  Flow  for  given  flight  parameters,  (5)  Water/Wind  Tunnel  Testings 
for  solution  validation,  and  (6)  Perform  3 DOF  Simulink  with  POD/RSM  data  for  MAV 
Stability/Control  demonstration. 

In  order  to  perform  such  an  aerodynamic  development  and  control  demonstration  for  a  mMAV, 
we  select  the  FLRW  designed  by  the  University  of  Florida  and  shown  in  Figure  30  as  a  testbed. 
The  FLRW  has  a  reflex  chamber  and  has  a  size  of  6  inches  with  a  weight  of  52  grams.  It  carries 
a  color  video  camera  and  is  powered  by  a  Maxon  electric  motor  to  achieve  speeds  of  12  to  30 
mph.  Control  is  accomplished  using  two  independently  controlled  elevons  that  are  actuated 
symmetrically  and  anti-symmetrically  using  small  rotary  servos.  In  addition,  because  of  the 
flexible  structure  that  provides  a  passive  adaptive  washout  mechanism,  the  FLRW  can  change  in 
wing  camber  and  wing  twist  in  gusty  wind  conditions.  This  is  to  reduce  the  sensitivity  of  the 
FLRW  to  disturbance.  Because  of  the  flexible  structure  and  the  reflex  camber  which  could 
increase  the  complexity  of  the  flow  structures,  the  FLRW  is  a  challenging  MAV  design  for 
accurate  aerodynamic  prediction.  The  success  of  the  tasks  proposed  will  certainly  demonstrate 
that  the  enhanced-aerodynamic  methodology  is  a  powerful  aerodynamic  tool  for  the  selected 
FLRW  with  a  view  that  it  can  achieve  the  Phase  II  technical  objectives  of  section  8.1  for  a  more 
general  class  of  MAV. 
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8.2.1  Rigid  mMA  V  wing  aerodynamics 

•  Generalize  the  2D  work  developed  in  Phase  I  to  3D  flow  for  a  rigid  mMAV  wing. 

•  The  method  will  use  RANS/CFL3D,  and  we  will  adopt  the  same  transition  flow  model 
technique  as  that  created  for  the  2D  work. 

•  The  obtained  CFD  solutions  in  terms  of  transition  point  locations,  transition  lengths, 
surface  pressures,  effects  of  vortex  stretching  and  flow  Reynolds  stresses  will  be  verified 
with  experimental  results  in  8.2.4. 

8.2.2  Flexible  mMA  V  wing  aeroelasticity 

This  is  a  grand  task  on  a  brand  new  look  on  the  aeroelasticity  of  mMAV  wing.  The  basic 
aeroelastic  procedure  of  a  tightly  coupled  aerodynamic  CFD  and  structural  interaction  for  the 
mMAV  wing  has  been  thoroughly  worked  out  previously  by  Shyy  et  al  [33].  However,  the 
procedure  is  based  on  a  Navier-Stokes/laminar  flow  model  without  the  consideration  of  flow 
transition.  Our  proposed  tasks  are  to  reconstruct  Shyy’s  procedure  using  the  proposed  unsteady 
RANS,  or  URANS/CFL3D,  and  to  further  generalize  the  rigid  aerodynamics  including  the 
transition  flow  technique  developed  in  8.2.1  to  that  for  a  flexible  mMAV. 

•  Reconstruct  Shyy’s  procedure  but  using  URANS/CFL3D  instead 

Moving  grid  technique 
Membrane  structural  solver 

•  Implement  the  transition  flow  model  technique  of  8.2. 1 

•  Select  computed  cases  for  comparison  (e.g.,  in  terms  of  rigid  to  elastic  ratio) 

Shyy’s  result  vs  present  result  (all  aeroealstic) 

Present  rigid  solution  versus  present  aeroelastic 

8.2.3  Gust  Consideration 

Adding  the  gust  model  into  the  rigid  and  flexible  mMAV  wing  is  an  essential  step  for  the  present 
3D  development. 

•  Validate  the  formulated  computational  Gust  approach/solutions  with  that  of  the  classical 
aeroelasticity,  i.e.,  Sharp-edge  gust  (Kussner  function)  and  Traveling  delta  function  gust 
(Sears  function),  etc. 

•  Generalize  the  2D  gust  handling  procedure  for  a  3D  rigid  mMAV  wing  (of  8.2.1). 
Evaluate  the  effect  on  flow  transition  and  separation  due  to  gust. 

•  Extend  the  gust  handling  procedure  of  8.2.1  to  a  flex  mMAV  wing  (of  8.2.2).  Evaluate 
the  flow  transition  and  flow  separation  effects  due  to  gust. 

•  Select  computed  cases  for  comparison  (e.g.,  in  terms  of  rigid  to  elastic  ratio). 

8.2.4  Wind/Water  -  Tunnel  Testing 

Professor  Luis  Bernal  of  the  ZONA  team  will  use  the  water-tunnel  and  the  wind-tunnel  facilities 
at  University  of  Michigan  for  testing  the  rigid/flexible  mMAVs.  Anticipated  steps  include 

•  Test  model  fabrication:  Rigid  model,  flexible  model. 

•  Equipments/DPI V  procedure  set  up. 


38 


•  Investigation  on  flow  about  a  pitching  3D  rigid  MAV  wing.  Measure  unsteady  effects  on 
rigid-wing  airload/response. 

•  Investigation  on  flow  about  Flexible  membrane  MAV  wing. 

•  Measured  data  is  to  be  validated  by  ZONA  with  ZONA  computed  solutions,  e.g.,  in  terms 
of  forces  and  moments. 


8.2.5  Perform  POD/RSM  of  the  Unsteady  Flow  of  3.2  for  Given  Flight  Parameters 

The  POD/RSM  development  here  is  to  provide  interpolated  solutions  for  stability  and 

aerodynamic  performance  evaluation  in  an  expedient  manner. 

•  POD/RSM  of  a  rigid  mMAV  wing  (of  8.2.1)  for  baseline  stability/flight  dynamic 
analysis. 

•  POD/RSM  of  a  flexible  mMAV  wing  (of  8.2.2)  for  stability/control  and  assure 
aerodynamic  performance. 

•  POD/RSM  should  treat  the  gust  modeling  as  an  individual  parametric  vector. 


8.2. 6  MA  V  Flight  Simulation  Model 

The  simulation  model  will  be  based  on  the  six  degree-of- freedom  nonlinear  body- fixed  equations 
of  motion  that  result  from  the  MAV  longitudinal  and  lateral/directional  decoupled  dynamics, 

.  Fx 

u  =  -  qw  +  rv  -  g  sin  6 

m 


Fv 

v  =  ^--ru  +  pw  +  geos  6  sin  cp 
m 

F 

w  =  -2-  —  pv  +  qu  +  g  cos  6  cos  (p 


m 

1  r 

P=T 


q=- 


L + Hxz  ~cir{fz  ~Iy)+ pq 4 

M-pr(lx-Iz)  +  r2Ixz-p2Ix 


q  = 


4 L 


N  ~  pq  ( 1  y  ~  Ix )  +  plxz  ~  qrlx 


(12) 


where  u(t),  v(t),  w(t),  p(t),  q(t)  and  r(t)  are  the  perturbed  forward,  lateral  and  vertical  speeds,  as 
well  as  the  roll,  pitch  and  yaw  rates  of  the  MAV,  respectively.  There  is  coupling  between  the 
longitudinal  and  lateral/directional  dynamics  due  to  the  inertial  and  gravitational  terms. 


The  MAV  physical  constants  to  be  considered  are  the  mass,  m,  the  wing  area,  S,  the  wing  span, 
b,  the  mean  aerodynamic  chord  c  ,  and  the  matrix  of  inertia  about  the  rigid-body  axis,  /,  that  is 


Ib 


(13) 
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The  aerodynamic  forces,  and  pitching  moment  are  produced  by  the  MAV  relative  motion  with 
respect  to  the  air  flow,  and  are  proportional  to  the  dynamic  pressure  qw  =0.5 pairV  as  well  as  to 
the  angle  of  attack,  a(t ) ,  and  sideslip  angle,  /3(t),  with  respect  to  the  relative  wind.  In  addition  to 
the  two  aerodynamic  angles  and  dynamic  pressure,  they  are  dependent  on  the  Re. 

Then,  the  aerodynamic  forces  (Fx,  Fy  and  Fzj  and  moments  (L,  M  and  N),  that  appears  in  the 
above  equation  will  be  computed  as  the  end  product  of  the  CFD  simulation  runs  by  the 
POD/RSM  approach.  Their  nonlinear  functional  dependence  on  the  aerodynamic  angles, 
dynamic  pressure  and  Re  will  be  modeled  along  this  Phase  II  effort  throughout  the  use  of  “look¬ 
up”  tables  in  the  Matlab/Simulink  software  environment. 


Figure  34  MAV  6-DOF  Implemented  in  the  Simulink  Environment. 

The  nonlinear  equations  of  motion  and  the  aerodynamic  forces  and  moments  computed  using  the 
proposed  POD/RSM  approach  will  be  implemented  in  the  MATLAB/Simulink  environment 
system  as  shown  in  Figure  34. 

It  is  clearly  observed  that,  the  inputs  are  the  symmetrical,  de,  and  anti-symmetrical,  5au,  elevon 
deflection  angles,  as  well  as  the  rudder  deflection  angle,  5rud ,  while  the  outputs  are  the 
longitudinal  and  lateral/directional  perturbed  variables. 

8.2. 7  Planned  Program  Schedule 

The  planned  program  schedule  for  performing  the  proposed  6  tasks  and  showing  the  distribution 
of  tasks  between  ZONA  and  University  of  Michigan  is  presented  as  follows: 
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Tasks 


Year  1 

Year  2 

1  2  3  4  5  6  7  8  9  10  11  12 

1  2  3  4  5  6  7  8  9  10  11  12 

To  Be 
Preformed 
By 


Task  l:Rigid  MAV  Wing  Aerodynamics 

•  Use  RANS  to  investigate  2D/3D  effects 

•  Use  RANS  to  obtain  Forces/Moments 

•  Validation  w/  Test  data 

Task  2:  Flexible  mMAV  Wing  Aeroelasticity 

•  Use  URANS  for  Cp,  Airload,  Response  on  wing 

•  Study  Shyy’s  results  vs  RANS  results 

•  Evaluate  Elastic/Rigid  Ratio 

•  Validation  w/  Test  data 

Task  3:  Gust  Considerations 

•  2D  Gust  models  revisit  and  CFD  formulation 

•  Use  URANS  for  2D  gust  effects  on  spt’n/transit’n 

•  Use  URANS  for  3D  gust  effects  on  spt’n/transit’n 

Task  4:  Wind/Water  Tunnel  Testing 

•  Equipment/  PIV  setup 

•  Investigate  the  steady  flow  of  a  3D  rigid  MAV  wing 

•  Investigate  flow  about  a  pitching  3D  rigid  MAV  wing 

•  Investigate  flow  about  Flexible  mMAV  wing 

Task  5:  Perform  POD/RSM  on  Data  of  Tasks  1/2 

Task  6:  Perform  6DOF/Simulink  for  mMAV  using 
Data  from  Task  5 

Task  7:  Final  Report  Documentation. 

Program  Management 

•Kick-off 
•Interim  Report 

•Final  Report/Final  Presentation 


A 


i  i  i 


ZONA/UM 


ZONA/UM 


ZONA 


UM 


ZONA 

ZONA 

ZONA/UM 


A 


Note:  UM  is  the  University  of  Michigan 
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LIST  OF  ACRONYMS  AND  ABBREVIATIONS 


ACRONYM 

DESCRIPTION 

AFRL 

Air  Force  Research  Laboratory 

ANSYS  ICEMCFD 

Commercial  CFD  meshing  program 

BSL 

Menter’s  Baseline  shear  layer  model 

CFD 

Computational  fluid  dynamics 

CFL3D 

NASA  CFD  solver  program 

Cp 

Pressure  coefficient 

DNS 

Direct  numerical  simulation 

FLRW 

Flexible  latex  rubber  wing 

J-L 

Jones-Launder 

LES 

Large  eddy  simulation 

Low-Re 

Low  Reynolds  number 

MAV 

Micro  air  vehicle 

mMAV 

Membrane  MAV 

POD 

Proper  orthogonal  decomposition 

R&D 

Research  and  development 

RANS 

Reynolds-averaged  Navier-Stokes 

RSM 

Response  surface  methodology 

S-A 

Spalart-Allmaras 

TU-BS 

Technische  Universitat  zu  Braunschweig 

UF 

University  of  Florida 

URANS 

Unsteady  RANS 

44 


